### A Pluto.jl notebook ###
# v0.20.5

using Markdown
using InteractiveUtils

# ╔═╡ f18a9488-13b8-11f0-342b-f99afaf3ac19
begin
	using Pkg
	Pkg.develop(path="/home/fritzsch/fri/JAC.jl")
end

# ╔═╡ fd4ef417-311f-43aa-a601-bf8a9d66b890
using JAC

# ╔═╡ 79374db5-1fbd-4499-b8d9-c7b840bdee52
html"""
<style>
	main {
		margin: 0 auto;
		max-width: 2000px;
    	padding-left: max(150px, 10%);
    	padding-right: max(150px, 10%);
	}
</style>
"""

# ╔═╡ 8d491602-2e68-4ce9-b949-859f91cfee9d
md"""
# Hydrogenic computations
"""

# ╔═╡ 75061d53-4fb7-45a6-b8bb-c199d52479da
md"""
Let us first invoke JAC, either from the own source-code basis or simply by **using JAC**
"""

# ╔═╡ 956d7bd6-d59f-4daa-88b6-9b8879918d79
md"""
Perhaps, the simplest atomic computations can be made for hydrogenic ions. If we consider, for example, hydrogen-like argon (Z=18), we can first determine and compare the energies of the ``1s`` and ``2p`` levels from the (non-relativistic) Schroedinger equation with those from Dirac's relativistic equation by
"""

# ╔═╡ b0f8586c-db7e-4872-b4e3-ff6beafb4c6c
begin
	Z     = 18.0
	e1s   = HydrogenicIon.energy(Shell("1s"), Z)   
	e2p   = HydrogenicIon.energy(Shell("2p"), Z)
	e2p_1 = HydrogenicIon.energy(Subshell("2p_1/2"), Z)   
	e2p_3 = HydrogenicIon.energy(Subshell("2p_3/2"), Z)
end

# ╔═╡ 4c2b4d2a-4b34-4f6a-9690-6efc9f2d3fea
md"""
While the (one-electron) energies are displayed by the function HydrogenicIon.energy() in the default units (presently eV here and as could be overwritten by the user), all **computations are internally performed and returned always in atomic units.** This applies to all computations if not stated otherwise in the description of some particular function; indeed, the user-defined units mainly in the tabulation of results but are not returned. Most of these tabulations are generated by display method which print some table to screen but return nothing otherwise.

This clear distinction between the use of **atomic** and **(user-specified) default units:: can be seen easily from the output of the variables `e2p, e2p_1, e2p_3` above, although these units can be converted also quite easily by:
"""

# ╔═╡ 759da271-0adf-4ae7-9642-8620c4c34690
begin
	(e2p, e2p_1, e2p_3) 
	e1s_eV = convertUnits("energy: from atomic to eV", e1s)
	e1s_au = convertUnits("energy: from eV to atomic", e1s_eV)
	(e1s, e1s_eV, e1s_au)
end

# ╔═╡ 8e2bf169-f65c-4f8c-b2a2-8ec18dd1e35f
md"""
... and similarly also for other energy units as well as for other physical entities; search for `convertUnits()` in the *LiveDocs*.

From these energies, we can easily compute the fine-structure splitting of the ``2p`` level into the ``2p_{1/2}`` and ``2p_{3/2}`` (fine-structure) levels for hydrogen-like argon from above:
"""

# ╔═╡ 6caf8fa3-4b0f-4c75-b968-1890ab258170
begin
	e13 = e2p_1 - e2p_3
	e13_eV = convertUnits("energy: from atomic to eV", e13)
end

# ╔═╡ a17abbaa-8f5c-49c1-9269-4f21018e6510
md"""
Apart from the (single-electron) energies, we can generate also the radial orbitals, i.e. the ``P(r)`` in the non-relativistic theory or ``P(r)`` and ``Q(r)`` in the relativistic theory. Both, the relativistic and non-relativistic orbitals, can be calculated by using the general (and analytically well-known) solutions of either the Dirac's or Schrödinger's equation.

However, since all (radial) orbital functions are represented on some `Radial.Grid`, special care has to be taken in order to define a proper grid, and which affects also all subsequent numerical computations as well as the accuracy that can be obtained. To better understand the definition of the grid, search for `Grid`, which shows the internal definition and the constructors to define such a grid. In practice, there are three (logical) levels that are considered and realized for each grid:

(i) The phenomenological grid definition that is based on our physics understanding and intuition of which grid is appropriate to decribe a particular atomic property or process. Two currently implemented grid refers to an exponential grid (with exponentially increasing stepsize) as appropriate for many bound-state computations and a linear-logarithmic grid that start exponentially but becomes linear in its stepsize for large values of ``r``. Such a linear-logarithmic grid is typically needed to describe electron ionization and capture processes as the incoming or outgoing electron waves are sinusoidal and should be represented with a proper number of grid points, say 20-30, per period of the electron wave. In JAC, the phenomenologically part of the grid is chosen (analogue as in GRASP) by the parameters `rnt`, `h`, `hp` and `NoPoints` in the definition above.

(ii) The phenomenological definition of the grid is internally translated into a **sequence of knots** upon which the B-spline primitives are defined. Here, each `nth` mesh point from the phenomenological is chosen to keep the number of B-splines and the corresponding size of matrices (that need to be diagonalized) moderate. These knots are used to define the B-splines and to determine the eigenvectors of all orbitals (within some given potential) but they are not applied in order evaluate matrix elements or radial integrals.

(iii) The phenomenological grid and the definition of the B-splines (knots) are eventually combined into a **physical grid** upon which all radial functions are represented. This physical grid resembles the phenomenological grid but with modified grid points in between the knots of the B-spline grid (t-grid). Here the radial points and corresponding weights are chosen due to a Gauss-Legendre distribution and goal to determine all radial integrals exact up to a given *Gauss-Legendre order*.

The physical grid is thus defined by the three arrays `r` (the mesh points along), `rp` (the derivatives ``d r/ dr`), `rpor` (the values `rp / r) as well as `w` (the corresponding weights). All these arrays are of length `NoPoints` but not necessarely equal as the grid points `nr` are coupled to the underlying integration, respectively, interpolation scheme.

In principle, this physical grid could chosen also on other interpolation/integration formulas, such as Gauss-Laguerre or others, that have been utilized in atomic physics. All what is needed would be to adapt the intermediate grid point and weights accordingly.

The clear distinction between the phenomenological grid and the physical grid help avoid that every new (radial) operator as well as particular boundary condition of the radial orbital functions need to be treated independently within the B-spline basis. Instead, a proper interpolation/integration formula should guarantee that all results are integrated sufficiently accurate, and this is first of all tested by enlarging the number of grid point `NoPoint --> nr`.

Here, we first apply an exponential grid, and which is appropriate below in order to compute various expectation values:
"""

# ╔═╡ 373afcb8-de1b-4aab-b14d-54a5f0db2fae
grid   = Radial.Grid(true, printout= true)

# ╔═╡ 3375b01f-3756-4ebf-9571-f3fcd0cf61db
md"""
With this grid, the non-relativistic radial ``P(r)`` orbitals can be obtained either for a single ``r``-value, for a list of ``r``-values as well as for all ``r``-values on a given radial grid:
"""

# ╔═╡ 2d66a365-5f6f-42c6-a97e-234f1fabc291
begin
	Pnr_1s = HydrogenicIon.radialOrbital(Shell("1s"), Z, grid)
	Pnr_2p = HydrogenicIon.radialOrbital(Shell("2p"), Z, grid)
end

# ╔═╡ 4573673b-6b51-46f9-bad9-f3c633bb3acc
md"""
Here, the exponential tails of the radial orbitals are simply set to zero if ``P(r) <`` `1.0e-15` (and similarly for ``Q(r)``).

Of course, we could plot the two functions `Pnr_1s` and `Pnr_2p` directly by some proper call of Plots, and if we give the correct arrays of radial mesh points from above (not shown yet).
"""

# ╔═╡ f371c925-a236-4962-bbd7-fdea676f2bf7
md"""
### Expectation values of hydrogenic ions
"""

# ╔═╡ 37e9cce9-7633-4fff-9b89-81abaee1c840
md"""
We have seen the radial dependence of the wave functions but shoud still check whether the generated wave functions are normalised, i.e. ``<1s|1s> \:=\: 1``. To do that, we can use the "Trapz" integration package, which simply performs trapezoidal integration over common Julia arrays. If you do not have this package installed, you can do this by running the following calls. Afterwards, the normalisation can be calculated simply by integrating the square of the wave function on the grid it was generated on. Note that for multiplying vectors, one needs to use ".*" instead of the simple multiplication sign.
"""

# ╔═╡ d1d38a33-44ad-449a-b49a-757f18b90281
begin
	Z1 = 1.
	Pnr_H_1s = HydrogenicIon.radialOrbital(Shell("1s"), Z1, grid)       
end

# ╔═╡ 8a5fbba2-c845-4077-b925-7f9a3c1a4e34
begin
	Pkg.add("Trapz")
	using Trapz
	norm = trapz(grid.r, Pnr_H_1s.*Pnr_H_1s)
end

# ╔═╡ a3bcfa37-1556-494c-b7f6-7080e933df99
md"""
We see that the normalisation is equal to unity, i.e. ``<1s|1s> \:=\: 1``, as expected. Further, we can verify other known properties of the hydrogen wave function. For example, the expectation values of the ``1/r,\; r`` and ``r^2`` operators are analytically known. As it is often taught in many university atomic physics courses, the expectation values are ``<1s| 1 / r |1s> \:=\: 1, a_o,\; <1s| r |1s> \:=\: 3/2\, a_o,\; <1s| r^2 |1s> \:=\: 3\,a_o``, where ``a_o`` is the Bohr radius. Let us verify this with JAC.
"""

# ╔═╡ 91abaa0f-3ec9-4b0a-9bca-4209c7cb1426
begin
	print(trapz(grid.r, Pnr_H_1s.*Pnr_H_1s./grid.r), "\n")
	print(trapz(grid.r, Pnr_H_1s.*Pnr_H_1s.*grid.r), "\n")
	print(trapz(grid.r, Pnr_H_1s.*Pnr_H_1s.*grid.r.*grid.r), "\n")
end

# ╔═╡ 1997bd27-c8fd-4cd8-8194-e5fad5767fa9
md"""
In JAC, all computations are internally performed and returned always in atomic units and hence ``a_o``. We can, therefore, see that our calculation confirms the analytically predicted expectation values.

Until now, we considered the solutions of the nonrealtivistic ion. In many cases, these nonrelativistic predictions are rather accurate. For heavy elements, however, the nonrelativistic theory comes short and the relativistic Dirac equation needs to be used instead. In the following examples, we will compare the nonrelativistic and relativistic radial wave functions and energies in order to get a intuition about the limits of the Schrödinger equation solutions. Firstly, let us generate the relativistic electron wavefunction for the hydrogen atom in its ground state, as well as the nonrelativistic and relativistic wave functions for hydrogen-like uranium.
"""

# ╔═╡ 854c4844-2ceb-451c-83ef-326587d30dbf


# ╔═╡ a9e25b61-ccd5-453c-bf91-a66e28de55f0


# ╔═╡ d70082f6-fa88-4b54-8306-46efe25138fc
md"""
# Still under construction !!!
"""

# ╔═╡ Cell order:
# ╟─79374db5-1fbd-4499-b8d9-c7b840bdee52
# ╟─8d491602-2e68-4ce9-b949-859f91cfee9d
# ╟─75061d53-4fb7-45a6-b8bb-c199d52479da
# ╟─f18a9488-13b8-11f0-342b-f99afaf3ac19
# ╠═fd4ef417-311f-43aa-a601-bf8a9d66b890
# ╟─956d7bd6-d59f-4daa-88b6-9b8879918d79
# ╠═b0f8586c-db7e-4872-b4e3-ff6beafb4c6c
# ╟─4c2b4d2a-4b34-4f6a-9690-6efc9f2d3fea
# ╠═759da271-0adf-4ae7-9642-8620c4c34690
# ╟─8e2bf169-f65c-4f8c-b2a2-8ec18dd1e35f
# ╠═6caf8fa3-4b0f-4c75-b968-1890ab258170
# ╟─a17abbaa-8f5c-49c1-9269-4f21018e6510
# ╠═373afcb8-de1b-4aab-b14d-54a5f0db2fae
# ╟─3375b01f-3756-4ebf-9571-f3fcd0cf61db
# ╠═2d66a365-5f6f-42c6-a97e-234f1fabc291
# ╟─4573673b-6b51-46f9-bad9-f3c633bb3acc
# ╟─f371c925-a236-4962-bbd7-fdea676f2bf7
# ╟─37e9cce9-7633-4fff-9b89-81abaee1c840
# ╠═d1d38a33-44ad-449a-b49a-757f18b90281
# ╠═8a5fbba2-c845-4077-b925-7f9a3c1a4e34
# ╟─a3bcfa37-1556-494c-b7f6-7080e933df99
# ╠═91abaa0f-3ec9-4b0a-9bca-4209c7cb1426
# ╟─1997bd27-c8fd-4cd8-8194-e5fad5767fa9
# ╠═854c4844-2ceb-451c-83ef-326587d30dbf
# ╠═a9e25b61-ccd5-453c-bf91-a66e28de55f0
# ╟─d70082f6-fa88-4b54-8306-46efe25138fc
